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ABSTRACT 

We  consider  the  solution  of  the  linear  systems  of  cJgebraic  equations  which 
cirise  from  elliptic  finite  element  problems  defined  on  composite  meshes.  Such 
problems  can  systematically  be  built  up  by  introducing  a  bzisic  finite  element 
approximation  on  the  entire  region  and  then  repeatedly  selecting  subregions,  cind 
subregions  of  subregions,  where  the  finite  element  model  is  further  refined  in 
order  to  gain  higher  accuracy.  We  consider  conjugate  gradient  algorithms,  and 
other  acceleration  procedures,  where,  in  each  iteration,  problems  representing 
finite  element  models  on  the  original  region  and  the  subregions  prior  to  further 
refinement  are  solved.  We  can  therefore  use  solvers  for  problems  with  uniform  or 
relatively  uniform  mesh  sizes,  while  the  composite  mesh  can  be  strongly  graded. 

In  this  contribution  to  the  theory,  we  report  on  new  results  recently  obtained 
in  joint  work  with  Mciksymili£m  Dryja.  We  use  a  basic  mathematical  frame  work 
recently  introduced  in  a  study  of  a  variant  of  Schwarz'  alternating  algorithm-  We 
establish  that  several  fast  methods  Ccin  be  devised  which  are  optimal  in  the  sense 
that  the  number  of  iterations  required  to  reach  a  certain  tolerance  is  independent 
of  the  mesh  size  as  well  as  the  number  of  refinement  levels.  This  work  is  also  tech- 
nically quite  closely  related  to  previous  work  on  iterative  substructuring  methods, 
which  cire  domain  decomposition  algorithms  using  non-overlapping  subregions. 

1.  Introduction.  In  this  paper,  we  consider  the  solution  of  the  large  linear  systems 
of  algebredc  equations  which  arise  when  working  with  elliptic  finite  element  approximations 
on  composite  meshes.  In  this  contribution  to  the  theory,  we  report  on  new  results  recently 
obtcdned  in  joint  work  with  MeiksymilieLn  Dryja. 

Finite  element  models  on  composite  meshes  cem  systematically  be  built  up  inside  a 
frame  work  of  conforming  finite  elements;  cf.  Ciarlet  [3].  We  do  so  to  be  able  to  use 
a  number  of  technical  tools  which  aie  available  primarily  in  the  conforming  case.  We 
begin  by  introducing  a  basic  finite  element  approximation  on  the  entire  region  and  we  then 
repeatedly  select  subregions,  cuid  subregions  of  subregions,  where  the  finite  element  model 
is  further  refined.  We  solve  the  resulting  linear  system  by  using  an  iterative  method  such  as 
the  conjugate  gradient,  Richardson  or  Chebyshev  method.  We  accelerate  the  convergence 
by  solving  so  cailled  standard  problems.  These  correspond  to  the  finite  element  models  on 
the  original  region,  prior  to  any  refinement,  Jind  those  on  the  subregions  prior  to  further 


*  CouTcint  Institute  of  Mathematical  Sciences,  251  Mercer  Street,  New  York,  N.Y. 
10012.  This  work  was  supported  in  part  by  the  National  Science  Foundation  imder  Grant 
NSF-CCR-8703768  and,  in  part,  by  the  U.  S.  Department  of  Energy  under  contract  DE- 
AC02-76ER03077-V  at  the  Courant  Mathematics  juid  Computing  Laboratory.  This  report 
was  prepared  for  the  proceedings  of  the  Second  International  Symposium  on  Domain  De- 
composition, which  wcis  held  at  UCLA  January  14-16,  1988. 


refinement.  We  ccin  thus  use  solvers  for  problems  with  uniform  or  relatively  imifonn  mesh 
sizes,  while  the  composite  mesh  Ccin  be  strongly  graded. 

This  approach  offers  a  number  of  advantages.  An  existing  code  can  be  upgraded,  in 
order  to  increiise  the  accuracy  locally,  without  a  radical  redesign  of  the  data  structures  etc., 
since  we  can  use  the  old  code  to  solve  one  or  several  of  the  standard  problems.  Issues  of  data 
structttres  and  geometry  are  generally  simpler  if  we  design  programs  for  composite  mesh 
problems  in  terms  of  simpler  steindaird  problems.  The  use  of  simple  standard  problems  also 
tends  to  improve  the  performance  of  the  programs  on  vector  machines.  Finally  we  note 
that  if  each  of  the  standard  problems  has  appreciatively  fewer  degrees  of  freedom  than  the 
composite  model,  then  we  might  benefit  from  solving  a  number  of  smaller  problems  rather 
than  one  large  one.  In  other  words,  we  might  view  our  approach  as  a  divide- and- conquer 
strategy. 

We  study  two  families  of  edgorithms  which  we  call  multiplicative  and  additive  respec- 
tively. The  so  called  additive  algorithms  are  pjirticularly  well  suited  for  parallel  computing 
in  that  in  each  iteration  step  we  can  simultaneously  solve  all  the  standard  problems  on  the 
the  different  levels  of  refinement.  Synchronization  between  the  processors  is  required  once 
in  each  iteration,  Ucimely  when  the  residual  is  computed  cind  assembled.  One  or  a  group 
of  processors  can  therefore  be  assigned  in  a  straightforward  way  to  each  of  the  standard 
problems.  As  we  wiU  see,  we  can  view  the  multiplicative  2dgorithin8  as,  preconditioned  con- 
jugate gradient  methods,  while  the  additive  variant  involves  the  solution  of  a  transformed 
equation  with  the  same  solution  &s  the  original  one  using  no  farther  preconditioning. 

The  systematic  study  of  the  methods  under  consideration  goes  back  at  least  to  1983, 
when  the  Fast  Adaptive  Composite  (FAC  )  method  was  introduced  by  McCormick  [12].  Is- 
sues related  to  the  implementation  of  this  method  on  parallel  computers  led  to  the  introduc- 
tion of  Asynchronous  FAC  (AFAC)  methods  [8]  a  few  years  later.  Numerical  experiments 
have  now  been  reported  for  model  cases  and  recently  Richard  Ewing,  Steve  McCormick  and 
others  have  begun  to  test  the  algorithms  for  more  diffictdt  problems  sirising  in  industry. 
The  convergence  of  the  FAC  method  is  discussed  in  McCormick  et  al.  [12],  [13]  under  a 
certain  additional  regularity  assumption.  An  important  contribution  to  the  theory  is  given 
by  Bramble,  Ewing,  Pasdak  and  Schatz  [2],  who  outlined  a  proof  of  optimality  for  a  two 
level  FAC  algorithm.  They  require  no  additional  regularity  in  their  proof.  In  recent  papers 
Bj0rstad  [1]  and  Maudel  and  McCormick  [11],  develop  a  theory  for  both  multiplicative  eind 
additive  algorithms,  using  two  levels.  In  pcirticular,  they  obtain  interesting  results  concern- 
ing the  relationship  between  the  spectra  of  the  two  methods.  In  a  recent  paper,  Mandel 
and  McCormick  [10]  establish  optimality  for  a  multi-level  AFAC  algorithm  for  a  special 
model  problem. 

Our  study  begim  with  the  discovery  that  FAC  and  AFAC  methods  have  a  structure 
quite  similar  to  that  of  the  classical  Schwarz  procedure,  see  Schwarz  [14],  aind  an  additive 
vsiriant  thereof  recently  considered  in  Dryja  [4]  and  Dryja  Jind  Widlund  [5].  Our  earlier 
work  was  in  turn  inspired  by  a  recent  paper  by  P.-L.  Lions  [9]  in  which  a  variational 
frame  work  for  the  classical,  multiplicative  Schwarz'  method  is  developed  for  continuous 
elliptic  problems.  We  were  able  to  establish  a  rate  of  convergence  which  is  independent  of 
the  number  of  degrees  of  freedom  as  well  as  the  ntmiber  of  subregions  for  a  special  kind 
of  additive  Schwarz  cdgorithm,  see  Dryja  [4]  and  Dryja  and  Widlund  [5].  In  oirr  view, 
the  central  theoretical  issue  of  the  iterative  refinement  algorithms,  which  are  discussed  in 
this  paper,  is  similar  to  that  of  Schwarz  type  algorithms  namely  the  design  and  study  of 
algorithms  for  which  the  rate  of  convergence  is  independent  of  the  nimiber  of  subproblems, 
i.e.  the  number  of  refinement  levels  as  weU  as  eind  the  mesh  sizes. 

In  section  2,  we  introduce  the  finite  element  problems  on  composite  meshes,  certain 


projections  and  ova  algorithms  in  their  basic  form. 

In  Section  3,  we  consider  a  multiplicative  algorithm  for  the  multi  level  case  and  show 
that  its  rate  of  convergence  is  optimal. 

In  Section  4,  we  introduce  several  multilevel  additive  algorithms  and  provide  a  nimiber 
of  bounds.  The  principal  result  is  that  one  of  these  methods  has  a  rate  of  convergence  which 
is  independent  of  the  nxmiber  of  refinement  levels,  eis  well  as  the  mesh  sizes. 

2.    Composite  Finite  Element  Problems  and  Basic  Iterative  Methods.   We 

consider  linezir,  self  adjoint,  elliptic  problems  discretized  by  finite  element  methods  on  a 
bounded  Lipschitz  region  ft  in  iZ".  To  simplify  the  presentation,  we  aiss\mie  that  the  differ- 
ential operator  is  the  Laplacian  and  that  we  \ise  continuous,  piecewise  linesir  finite  elements. 
However,  almost  all  of  our  results  can  be  extended  immediately  to  general  conforming  finite 
element  approximations  of  any  self  adjoint  elliptic  problem  which  can  be  formulated  as  a 
minimization  problem.  The  continuous  and  discrete  problems  are  of  the  form 

aiu,v)  =  f{v),^veV, 

cLnd 

a(uh,vh)  =  /(t;h),  Vt;h€F\  (1) 

respectively.  The  spaces  V  and  V^  are  defined  in  the  next  few  paragraphs.  The  bilinear 
form  is  defined  by 

a{u,  v)  =       Vu  •  Vv  dx  . 
Jn 

This  form  defines  a  semi-norm  \u\hi  =  (a(ii,  u))^/'  in  H^{Cl). 

To  simplify  the  presentation,  we  assimie  that  we  have  a  zero  Dirichlet  boundeiry  con- 
dition on  dfl,  the  boundary  of  f2.  The  space  V  is  thus  5^o(fi).  We  note  that  we  equally 
well  could  have  considered  an  inhomogeneous  Dirichlet  problem.  The  space  V  is  defined 
on  a  composite  triangulation,  which  is  possibly  the  resiilt  of  a  large  number  of  successive 
refinements.  The  triangulation  of  fi  is  given  in  the  following  way. 

We  first  introduce  a  relatively  coarse  tricingulation  of  fi,  eilso  denoted  by  fti,  amd 
denote  the  corresponding  space  of  finite  element  functions  by  V^^ .  We  can  think  of  this 
space  as  having  a  relatively  uniform  (or  uniform)  mesh  size  hi .  Let  CI2  be  an  area  where  we 
wish  to  increase  the  resolution.  We  do  so  by  subdividing  the  elements  auid  introducing  an 
additional  finite  element  space  V'*'.  We  cissure  that  the  resulting  composite  space  V''*  -I- 
V''*  is  conforming  by  having  the  functions  of  V^^  vanish  on  50?.  We  repeat  this  process  by 
selecting  a  subregion  CI3  of  1^2  ^^^  introducing  a  further  refinement  of  the  mesh  and  finite 
element  space  etc..  We  denote  the  the  resulting  nested  subregions  and  subspaces  by  n<  and 
V'^  respectively.  Throughout,  we  have  ^i  C  fii-i  and  V*^-^  n  H^itti)  C  V^  C  HliQi), 
i  =  2, . . . ,  Jb.  The  composite  finite  element  space  on  the  repeatedly  refined  mesh,  is 

V^  =  V'''  +V^^  +  ...  +  F''*  . 

We  assume  that  all  the  elements  are  shape  regular  in  the  sense  that  there  is  a  uniform 
bound  on  Hk/ Pk-  Here  hjc  and  px  are  the  diameter  and  the  radius  of  the  largest  inscribed 
sphere  of  any  element  K,  respectively.  Our  bounds  in  the  theory  developed  below  also 
depend  on  the  shape  of  the  subregions  Qi.  Thus  in  order  for  our  proofs  to  work,  we  cannot 
cdlow  the  sets  Qi_i  \  Cli  to  become  arbitrarily  thin  in  comparison  with  the  diameter  of  Cli-i. 
We  also  assume  that  the  £irea  of  any  triangle  on  level  i  can  be  bounded  by  const.  q^~*,  j  <  i 
times  the  area  of  the  triangle  on  level  j  of  which  it  is  a  part.  Here  9  is  a  constant  <  1. 


The  finite  element  problem  is  defined  by  equation  (1)  and  the  corresponding  stif&iess 
matrix  can  conveniently  be  computed  by  using  a  process  of  subassembly.  Introducing 
subscripts  to  indicate  the  domain  of  integration,  we  write 

a(u,  v)  =  an,\n,(",  v)  +  an,\n, (u,  v)  +  . . .  +  oo»  (tx,  v)  . 

The  stif&ess  matrices  corresponding  to  the  regions  Cli  \  fii+i,  i  <  k  —  1,  and  Clk  are 
computed  by  working  with  basis  functions  related  to  the  mesh  size  hi  .  The  quadratic 
form  corresponding  to  the  composite  stii&iess  matrix  is  the  stmi  of  the  quadratic  forms 
corresponding  to  fij  \  Qi+i  and  fife.  In  our  algorithms,  we  use  solvers  for  the  same  elliptic 
problem  on  the  subregions  Hi  ,  i  =  1, ...,  k,  and  the  relatively  uniform  meshes  corresponding 
to  hi  .  The  corresponding  stii&ess  matrices  can  in  fact  be  obtadned  at  a  small  extra  cost 
during  the  assembly  of  the  composite  mesh  model.  When  we  refine  a  finite  element  model 
locally,  the  modified  stii&ess  matrix  is  obtained  by  replacing  the  quadratic  form  associated 
with  the  subregion  in  question  by  the  one  corresponding  to  the  refrned  model  on  the  same 
subregion.  It  is  therefore  relatively  easy  to  design  a  method  which  systematically  generates 
the  stiffiiess  matrices  for  all  the  stiindard  problems  necessary  while  the  sti&ess  matrix  of 
the  composite  model  is  computed. 

The  fundamental  building  blocks  of  our  cilgorithms  are  the  PJ,  i  <  j,  the  projections 
onto  the  spaces  V^'  f] Hq{Q,j).  We  note  that  if  j  >  i,  then  we  solve  »  problem  on  fi^  with 
a  cojirser  mesh  than  if  V**'  were  used.  The  projection  PJ,  i  <  j,  is  defined,  in  terms  of  the 
tmique  element  of  F***  D  Hq{SIj),  which  satisfies 

a{P;vH,4>h)  =  a{vH,4>H),y4>heV''\  (2) 

We  now  introduce  the  multiplicative  and  additive  cilgorithms.  Since  no  further  effort 
is  involved,  we  develop  a  framework  which  is  also  useful  in  other  contexts  such  as  the 
study  of  algorithms  of  Schwaxz  tjrpe;  cf.  Dryja  [4],  Dryja  and  Widlund  [5]  and  Lions 
[9].  We  begin  by  considering  the  Ccise  of  two  subspaces  Vi  and  Vj  and  the  multiplicative 
(sequential)  algorithms.  In  a  first  fractional  step,  we  find  a  correction  ^lu"  G  Vi  of  the 
current  approximation  u"  by  solving 

aiSivT^v)  =  f{v)  -  a(u",t/)  =  a{u*  -  u",t»),  V  w  e  Vi. 

Here  «*   G  Vi  +  Vj  is  the  solution  of  the  given  problem.    The  calculation  of  a  second 
correction  ^ju"  €  Vj  completes  the  (n  +  l)th  step. 

a{62vr,v)  =  f{v)  -  a(u"  +  *iii",  v)  =  a{u*  -  («"  +  tfju"),  v),  Vr  6  Vi  . 

As  shown  in  Lions  [9],  it  is  easy  to  see  that 

6xvr  =  Pi[xi.' -  vT) 

SiU"  =  P2(u*  -  u"  -  5iu")  =  Piil  -  Pi)(«'  -  u") 

and  thus  the  error  propagates  as 

tx"+^  -u*  =  (J-P2)(7-Pi)(u'*-U*). 

Here  Pi  £ind  Ft  aie  the  orthogonal  projections  associated  with  the  bilineeir  form  a(-,  •)  and 
Vi  and  Vi,  respectively. 


We  can  thus  view  this  algorithm  as  a  simple  iterative  method  for  solving 

{Pi  +  P2  -  P2Pl)uH  =  9h, 

with  an  appropriate  right  hand  side  qh-  We  note  that  this  operator  is  a  polynomied  of  degree 
two  2ind  thus  not  ideal  for  parallel  computing  ,  since  two  sequential  steps  are  involved  . 
If  we  use  more  than  two  subspaces  and  therefore  more  projections  this  effect  is  further 
pronoimced.  The  basic  idea  behind  the  additive  form  of  the  algorithm  is  to  work  with 
a  simplest  possible  polynomial  in  the  projections.  With  k  subspaces  we  thus  solve  the 
equation 

Puh  =  iPi+P2  +  ...  +  Pk)uH  =  g'n  ,  (3) 

by  an  iterative  method.  If  we  can  show  that  the  operator  P  is  sjrmmetric  and  positive 
definite,  the  iterative  method  of  choice  is  the  conjugate  gradient  method  at  least  on  a 
computer  with  conventional  architecture.  We  must  also  make  sure  that  this  equation  has 
the  same  solution  Jis  equation  (1),  i.e.  we  must  find  the  correct  right  hand  side.  Since  by 
equation  (1),  we  have 

a{uh.,4>h)  =  f{4>h)  , 

we  can  construct  the  right-hand  side  5^  by  solving  this  equation  restricted  to  all  the  different 
subspaces  and  adding  the  results.  It  is  similarly  possible  to  apply  the  operator  P  of  equation 
(3)  to  any  given  element  of  V''  by  applying  each  projection  Pi  once.  Most  of  the  work  ,  in 
particular  that  which  involves  the  individual  projections,  can  be  carried  out  in  parallel. 

It  is  weU  known  that  the  number  of  steps  required  to  decrease  an  appropriate  norm  of 
the  error  of  a  conjugate  gradient  iteration  by  a  fixed  factor  is  proportional  to  ^/k,  where 
K  is  the  condition  ntmiber  of  P;  see  e.g.  Golub  Jind  Vaui  Loan  [7].  We  therefore  need  to 
establish  that  the  operator  P  of  equation  (3)  is  not  only  invertible  but  that  satisfactory 
upper  and  lower  bounds  on  its  eigenvalues  Ccin  be  obtained. 

We  end  this  section  by  discussing  a  symmetric  version  of  the  multiplicative  algorithm. 
The  projections  generally  do  not  commute  and  if  we  wish  to  have  a  sjrmmetric  expression  in 
the  operators  Pi,  which  allows  us  to  accelerate  the  convergence  by  the  stzmdard  conjugate 
gradient  method,  we  have  to  use  additional  fractional  steps.  In  the  case  of  two  subspaces, 
we  can  solve  the  first  problem  again.  Since  (/  -  Pj  )*  =  (/-  -P:)  we  can  write  the  resulting 
operator  as 

/  -  (/  -  Pi)(/  -  P2){I  -  Pi)  =  Pi  +  P2  -  P1P2  -  P2P1  +  P1P2P1  =1-  T2T;, 

where  T2  =  (J  -  Pi) (J  -  Pj).  The  corresponding  operator  for  ib  subspaces  involves  2^-1 
fractional  steps  and  has  the  form 

I-nn,  whereTfe  =  (J-Pi)(7-P2)---(J-Pfc). 

The  error  propagation  operator  for  the  basic  multiplicative  algorithm  is  T^,  while  the 
convergence  rate  of  the  symmetrized  algorithm  can  be  bounded  in  terms  of  the  condition 
number  of  J  -  TkTJ^.  We  note  that  a  bound  on  the  spectral  radius  of  Tfc  is  obtained 
inamediately  from  the  spectral  bounds  on  the  symmetric  operator. 

3.  The  Multiplicative  Algorithm.  In  this  section,  we  establish  the  following  restilt. 
We  note  that  our  aneilysis  resembles  that  of  Bramble  et  al.  [2]  in  the  case  of  fc  =  2. 

Theorem  1.  The  symmetrized,  multiplicative  iterative  refinement  algorithm  based  on 
the  projections  P/ ,  i  <  k,  has  a  condition  number  which  is  independent  of  k  and  the  number 


of  degrees  of  freedom.  The  spectral  radius  of  the  basic  multiplicative  algorithm  is  bounded 
by  a  constant  which  is  uniformly  less  than  one. 

We  begin  by  designing  a  preconditioner  b{u,v)  for  the  finite  element  problem  on  the 
space 

V'' =  V''' +'••  +  ¥'"'. 

We  then  show  that  this  quadratic  form  can  be  bounded  uniformly  from  above  and  below 
by  the  quadratic  fonn  of  the  composite  finite  element  problem.  Finally  we  establish  that 
we  can  work  practically  with  this  preconditioner  and  that  the  error  propagation  matrix  is 
the  same  as  that  of  the  S3rmmetrized  multiplicative  edgorithm,  introduced  in  the  previous 
section.  The  algorithms  are  thus  the  same. 

In  addition  to  the  projections  P?,  we  also  use  another  family  of  operators,  H'^,  i  = 
j  -  l,j,  defined  by 

Hivnix)  =  Vfc(x),    X  G  fi  \  fi,-, 
a{HJvH,  wh)  =  0,    V7i;h  €  V"^  n  H^iQj). 

We  can  call  HjVh  the  /ij— harmonic  extension  of  t;^  to  fij,  since  it  is  the  solution  in 
F''*  of  a  discrete  Dirichlet  problem  with  zero  right  hand  side  and  with  boundary  data  on 
ddj  given  by  vs- 

In  order  to  prepare  for  the  work  that  remains,  we  formulate  three  lemmas^ 

Lemma  1.  HiZiSi~^  =  HlZl,  3  <  i  <  *. 

The  proof  foUows  immediately  from  the  definition. 

Lemma  2.  a(5,'-^Uh,Fr^t;h)  =  a(irr^UH,Uh),  Vu^  e  V",  Vv^  G  V^*-\  2  <  i  <  k. 

Proof:  From  the  definition  of  HJ'^  foUows  that  Hi'^Vh-vn  G  V**-*  nH^{Qi).  Since 
Hl~^Uh  is  /i<_i— heirmonic  on  fij  the  result  follows  by  the  orthogonality  inherent  in  the 
definition  of  /ij_i— harmonic  functions. 

The  following  lemma  is  a  consequence  of  the  extension  theorem  given  in  Widlund  [15]. 

Lemma  3.  There  exists  a  constant  C  which  is  independent  of  u^  and  the  mesh  sizes, 
such  that  for  i  =  j  —  1,  j, 

an,{H}uh,H}uH)  <  Can,_.\n,(«H,«H),  V«h  G  V^. 

We  will  not  discuss  the  proof  of  this  lemma.  We  note  that  the  constant  C  necessarily 
blows  up  if  we  let  the  area  of  f2j_i  \  fij  shrink  to  zero,  keeping  fij-i  fixed.  Such  situations 
are  not  of  particular  interest  in  our  applications. 

Trivially,  we  can  write  Uh  ~  Pk^h  +  {I  -  -P* )"h'  I*  is  easy  to  see  that  the  second  term 
equals  HJ^Uh.  The  two  terms  are  orthogonal  in  the  sense  of  the  bilinear  form  and  thus 

a(uH,  vh)  =  a{PJ:uH,  Pfc^h)  +  a((/  -  P*K,  (/  -  P*K) 
=  a(Pfc*Uh,P*t;H)  +  a(5*Uh,ir^h)  . 

We  now  introduce  a  preconditioner,  which  in  the  case  k=2  is  the  final  one  but  which  in  the 
general  case  is  only  a  first  step  of  our  construction. 

aiPHuK^PllvH)  +  aiHJ:-'uH,HJ:-'vh)  . 


The  /ifc— hsumonic  function  on  fife  has  thus  been  replaced  by  the  /ife_i— harmonic  function 
with  the  same  bound£iry  values.  Since  the  latter  space  of  discrete  harmonic  functions  is 
smaller,  it  is  easy  to  see  that  the  preconditioner  is  bounded  from  below  by  the  original 
quadratic  form 

To  find  a  bound  from  above,  we  have  to  show  that  the  energy  attributable  to  the  hk-i- 
harmonic  function  can  be  bounded  by  that  of  the  /ifc— hcirmonic  function.  This  follows 
directly  from  Lemma  3. 

When  we  now  turn  to  the  case  of  /b  >  2,  we  repeatedly  replace  certeiin  discrete  harmonic 
functions  by  others.  In  each  step,  we  replace  the  current  last  term  of  the  preconditioner  by 
two.  We  begin  by  considering  the  identity 

H',-'u,  =  P',:lHl:-W  +  (/  -  PHlDHt'riH  . 

It  is  easy  to  see  that  the  second  term  equals  H^Zl ^k~^ ''^h-  As  in  the  first  step,  we  replace 
this  last  term  by  another,  namely  H^zl^k'^^'^i  *^*^  ^®  obtain  a  preconditioner  with  three 
terms: 

We  simplify  this  expression,  by  using  Lemma  1  2md  replace  the  third  term  of  the  precon- 
ditioner by  a{H^Zi''J^h, S^Zi'fh)-  By  repeating  the  process  just  outlined,  we  curive  at  the 
following  preconditioner: 

b{uH,VH)  =  a(P^H,Pfc^h)  +  a{PJ:zlHJ:-'u^,Pl:zlHJl-'vf,) 
+  • .  •  a{PlHluH,  PIhIvh)  +  a{H\uH,  H\vh). 

We  can  now  complete  our  proof  of  the  optimality  of  the  preconditioner.  By  using  the 
definitions  of  the  operators  P?  and  Hp  we  find  that 

aiP^UH^Ptun)  =  an,(Pfc*Uh,Pfe*Uh) 
and,  for  3  <  i  <  is, 

a{PiziHr'nH,pizlsr'u^  =  a^,_,{PiziHir'uH,piziHi-'uH) 
=  ao,_.(^r'«H»^r'«h)  -  an._AsiziHr'u^,H:ziHr'uH) 

=  ao,_ASr'uH,Hr'uH)  -  ao,.,{HtzluH,Hizl^H)^ 

In  the  last  step,  we  have  used  Lemma  1. 

By  using  the  definition  of  the  Hj  operators,  we  can  rewrite  the  preconditioner  as 

H'^h,Vh)  =  an»(uh,uh)  -  anJ5'fcUh,5'fcUh) 


This  quadratic  form  can  be  written  as  the  sum  of  an(tihj  Uh)  and  a  number  of  terms  of  the 
form 

Since,  on  fii,  H\~  Uh  and  Hluh  are  /ij_i— haumonic  and  /ij— harmonic  functions,  respec- 
tively, with  the  same  boundziry  values  on  ^fij,  it  is  eaisy  to  see  that  these  terms  are  positive. 
This  shows  that  the  preconditioner  is  bounded  from  below  by  an(u/,,«h).  To  get  an  upper 
bound,  we  only  have  to  estimate  the  positive  terms.  By  Lemma  3,  the  energy  attributable 
to  the  subregion  fij,  of  the  /ij_t— harmonic  function  El~^Uh  can  be  estimated  by  a  constant 
times  ani_i\nj(i^h»"h)'  The  proof  of  the  upper  botmd  is  completed  by  simiming  over  i. 

What  remains  is  to  establish  that  the  use  of  this  preconditioner  leads  to  a  series 
of  problems  which  are  directly  related  to  the  projections  P*  in  particular  the  S3rmmetric 
multiplicative  algorithm  given  in  section  2.  For  the  unaccelerated  algorithm,  the  equation 
that  we  have  to  solve  is  of  the  form 

h{8uh,VH)  =  a{eh.,VH)  , 

where  e^  is  the  error  before  the  current  step  and  Suh  th^  next  correction.  We  first  choose 
test  functions  v^  in  the  subspace  F'**.  All  terms  of  the  preconditioner,  except  the  first, 
vanish  for  such  test  fimctions.  Therefore  the  result  of  this  first  fractional  step  equals 
P^Suh  =  Pk^h'  We  choose  V^*-'  as  our  second  space  of  test  functions.  All  the  terms 
of  the  preconditioner  except  the  first  two  vanish  cind  the  first  is  alreeidy  known  and  can 
therefore  be  moved  over  to  the  right  hand  side.  A  straight  forward  Ccilailation  and  Lemma 
2  show  that  the  new  right  heind  side  is  equal  to  a((l  —  Pll)eh,Vh).  In  this  second  fractional 
step,  we  obtain  PJ^Zi  ^k"^ ^^^  =  PJ^Zli^  —  Pk)^h'  Proceeding  in  this  manner,  we  compute, 
in  the  *-th  fractional  step,  HlSuh  =  Pi{I  —  Pi)  •••{!—  Pk)eh-  At  this  time,  the  values  of 
Stth  are  available  on  fii  \  fij-  We  can  use  its  boundary  values  on  dfli  as  Dirichlet  values  and 
solve  a  problem  in  V'*'  to  obtain  the  values  of  Suh  in  fij  \  fia.  Proceeding  in  this  manner, 
using  a  total  of  (2A:  —  1)  fractional  steps,  we  finadly  obtain  Sun  everywhere.  A  calculation 
shows  that  the  error  propagation  operator  corresponding  to  the  whole  step  is  equal  to 

This  shows  that  the  method  defined  by  the  preconditioner  indeed  is  the  same  aa  that  of 
the  symmetrized  multiplicative  algorithm  defined  in  section  2. 

4.  Some  Additive  Algorithms.  We  recall,  that  in  the  so  called  additive  algorithms, 
we  solve  equation  (3)  by  the  conjugate  gradient  or  some  other  standard  iterative  algorithm. 
The  different  algorithms  can  simply  be  defined  by  specifying  the  subspaces  Vi  or  alterna- 
tively the  projections  Pi .  Before  we  discuss  specific  algorithms,  we  make  some  remarks  on 
estimating  the  eigenvalues  of  P  from  above  and  below. 

The  upper  bound  on  the  spectrum  is  obtained  by  bounding 

a{Pvh,Vh)  =  a{PiVh,VH)  +  a{P2VH,VH)  +  • '  •  +  a{PkVf„VH)  , 

from  above  in  terms  of  a{vh,Vh)-  We  can  use  Schwjirz'  inequality  and  the  fact  that  Pi  is 
a  projection  to  prove  that  each  term  is  bounded  by  a{vh,  Vh)  iind  thus  the  spectrum  of  P 
is  boimded  from  above  by  k.  Our  goal,  however,  is  to  establish  a  uniform  bound  on  the 
condition  number.  This  can  be  done  if  the  terms  are  orthogonal  or  almost  orthogonal;  cf. 
discussion  below. 


A  lemma  from  Lions  [9]  provides  a  method  for  obtidning  lower  boimds.  Since  the  proof 
of  his  result  is  quite  short,  we  include  it  in  this  paper. 

Lemma  4.  Lei  Uh  =  l^j_i  Uh,u  vihere  Uh,i  €  Vi,  be  a  partition  of  an  element  ofV*^  and 
assume  fuHher  that  ^i-^  a{uh,i,Uh,i)  <  Coa(«h,«h),  Vuh  e  V*.  Then  Xmin{P)  >  C^^. 

Proof:  By  elementairy  properties  of  s)rQimetric  projections  and  the  representation  of 
Uh  as  a  sum,  we  find  that 

k  k  k 

a{uh,Uh)  =  ^a{uh,Uh,i)  =  '^  a{uh,  PiUh,i)  =  ^  a(PjUh,Uh,i)  . 

t=l  »=1  i=l 

Therefore, 

fc  k 

a(uh,«h)  <  {'£a{PiUH,PiUH))'^'{J^aiuH,uUH,i)y^'  . 

t=i  t=i 

By  the  cissumption  of  the  lemma 

fc  fc 

a("h,«h)  <  Cl  ^  a{PiUH,PiUH)  =  Cl  ^  a(P.-Ufc,Uh)  =  Cla{PuH,Uh), 

and  the  lemma  is  established. 

We  consider  three  different  algorithms  eind  distinguish  between  them  by  using  a  super- 
script 1,  2  or  3.  The  most  naturzil  algorithm  cimoimts  to  iising  the  projections  f*,-  =  Pi, 
where  the  projection  operators  P,*  have  been  defined  in  equation  (2).  The  condition  number 
of  this  cilgorithm  grows  linearly  with  ik.  We  have  adready  remarked  that  the  eigenvalues  of  P 
always  are  boimded  from  above  by  k.  This  boimd  is  attained  if  V''*  PI  Hq  (fife)  is  not  empty, 
i.e.  when  the  mesh  size  is  fine  enough.  Any  such  function  belongs  to  V^\  i  =  1,2, . .  .,i, 
and  is  exactly  reproduced  by  each  of  the  projection  operators.  It  is  thus  an  eigenfunction 
with  eigenvalue  k.  Similarly,  any  function  which  belongs  to  V**'  n  Hq{^i  \  fij)  is  an  eigen- 
function with  eigenvalue  1.  We  will  show  later  that  all  the  eigenvalues  are  bounded  from 
below  by  a  constant,  which  completes  the  proof  that  the  condition  number  of  P^^'  is  of 
order  k. 

A  very  promising  method,  for  which  to  our  knowledge  the  optimality  has  only  been 
established  in  a  quite  specifd  model  case,  cf.    [10],  is  based  on  using  P,-      —  P,'  -  P,*^.i  , 

t  <  Jb  -  1  and  P^'^  =  P*  au  the  basic  projections  in  equation  (3).  It  easy  to  show  that  these 
differences  of  projections  axe  projections  and  that  the  composite  finite  element  space  V 
is  the  direct  sum  of  the  corresponding  subspaces.  A  difficulty  experienced  when  trying  to 
establish  that  the  eigenvalues  of  P^^^  are  bounded  uniformly  from  below,  by  using  Lions' 
lemma,  is  related  to  this  complete  lack  of  flexibility  in  representing  a  given  element  of 
y''  as  a  sum  of  elements  of  the  subspaces.  So  fax  we  have  only  been  able  to  prove  that 
Cq  <  const.  I. 

We  can,  however,  prove  the  following  result. 

Theorem  2.   The  eigenvalues  ofP^^^  are  uniformly  bounded  from  above  by  a  constant. 

Before  we  turn  to  the  proof  proper,  we  observe  that  the  unbalance  of  the  first  method, 
which  resulted  in  the  cimplification  of  certain  functions  by  a  factor  k,  is  no  longer  possible. 
By  regrouping  terms,  we  also  see  that 

p<^)  =  Pi  +{p^-p^)  +  ...  +  {p>;  - p;;-' ) .  (5) 
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All  the  terms,  except  the  first,  are  similar  to  mtiltigrid  corrections  since  they  each  represent 
the  difference  between  two  solutions  on  the  same  subregion,  using  two  different  mesh  sizes. 
By  using  the  representation  of  P^^\  given  in  equation  (5),  we  obtain 

p^^^uH  =  PluH  +  [p!  -p^)uh  +  ---  +  {pH  - p*-^ K 

=  Ui  +  U2  H +  Ufc. 

We  show  that  these  terms  aire  increasingly  orthogonal: 

a(u£,u„)  <  const.  gj'-'"l(a(ti,, ti,))^/*(a(ti„,u^))l/^ 

where  qi  <  I,  uniformly.  This  is  a  so-called  strengthened  Cauchy  inequality.  Without 
loss  of  generality,  we  assume  that  t  <  m.  It  is  easy  to  show  that  a(um,Vh)  =  0,  Vv^  £ 
yf^-i  n  H^{Qrn),  i-e.  u„  is  /i,„_i -harmonic  on  n„,.  Let  u<  =  Sa,</>i'  where  the  (?^|"  are 
basis  functions  in  V^' .  Any  such  bsisis  function  is  orthogonal  to  u„  if  its  support  does  not 
intersect  dClm  since  then  either  its  support  belongs  to  fi<  \  fl^  or  <t>'^'  £  V'^-^  n  HQ{Um)' 
Thus  the  only  contributions  to  a{um^u^)  originate  from  the  triangles  of  F'*'  which  intersect 
dQm-  Consider  one  such  triangle  T.  The  restriction  of  u/  to  T  is  a  linear  function  and  thus 
has  a  constcint  gradient.  It  can  also  be  written  as  a  linear  combination  of  basis  functions 
in  V''"'-'.  The  support  of  most  of  these  do  not  intersect  5n„,  and  therefore 

aT(«m,Ui)  =  arn5(«m,«<)  <  (ar(«m,«m))^^^(arn5(ti£,«<))^''^ 
<  Const.  gl'"~"ar(u„.„)'/'aT(u,,u^)^/^  . 

Here  S  denotes  the  union  of  the  small  trijuigles  that  are  next  to  dClm'  We  also  use  o\ir 
assumptions  on  the  triangulations,  given  in  section  2,  and  the  fact  that  the  gradient  of 
ut  is  constant  on  T  and  that  therefore  its  contribution  to  the  quadratic  form  is  directly 
proportional  to  the  eirea  of  integration.  The  proof  of  the  strengthened  Cauchy  inequality 
is  completed  by  summing  over  the  triangles  of  V**'. 
We  now  note  that 

k 

a{P^'^UH,P^'^UH)  =    J2  '^("''«^)  ^  ^^^^• 
i,3=l 


Here 


1  qi          qi  "•  \ 

A  =  const.  \     ^'          ^          «i  •••  «i"' 

i.            .            .  ...  . 

„k-l  „h-t         k-3  1 

9i  ?i  9i  •••  1 


cind 


It  is  easy  to  show  that  the  Euclidean  norm  of  A  is  imiformly  bounded.  Thus 

k  k 

a(p(2)u^p(2)„)  <  \A\t,Y^aiui,Ui)  <  const.  2a(u<,m). 

To  complete  the  proof,  we  note  that 

a{ui,Ui)  =  a{PluH  -  Pr'uH^PluH  -  Pf^H) 
=  a(P>h-P'-^Uh,P>h), 
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since  P'tx^  —  P,'     Uh  it  /ij_i— harmonic.  By  using  elementary  properties  of  the  projections, 
we  find  that  a(ui,ti,)  =  a{ui,Uh).  Therefore 


^  a{ui,  Ui)  =  ^  a{ui,  uh)  =  a{P^^Kh,  uh)- 


Thus 

a(p(2)uh,p(2>Uh)  <  const.  a{P<'^UH,UH), 

from  which  the  upper  bound  on  P^^^  follows. 

If  we  modify  our  projections,  we  airive  at  an  algorithm  for  which  we  have  a  proof  of 
optimality. 

Theorem  3.  The  additive  method  defined  by  the  projections  P,-  '  =  Pi~Pi+i,  i  <  *-2, 

P^_j  =  PJ^Zl  and  P^  —  Pj^  has  a  condition  number  which  is  independent  of  k  and  the 
number  of  degrees  of  freedom. 

A  uniform  upper  bound  is  obtcdned  by  using  Theorem  2  and  the  following  formula. 

p<')  =  Pi  -  p^  +PI-PI  +  ...  +  p*r/  +  p* 
=  pi-p^  +  pi-pi  +  '--  +  pl: 

+  p^-pl+p^-Pl  +  ...  +  pl;:l 

-  pW  +  p(2) 

By  Theorem  2,  there  is  an  upper  bound  for  P^^^  and  the  same  result  also  provides  a  uniform 
upper  bound  for  P^^\  which  is  an  operator  corresponding  to  a  composite  finite  element 
problem  using  k  —  I  levels.  ^ 

From  this  argument  juid  the  fact  that  P^^^  is  positive  definite,  it  follows  that  the 
condition  number  of  P^'^  cannot  be  much  better  than  that  of  P^^K  A  lower  botmd  of  the 
spectrum  of  P^^^  is  given  by  a  partitioning  formula  and  Lions'  Lemma.  Before  we  give  the 
details,  we  note  that  this  argument  also  provides  a  lower  bovmd  for  the  spectrum  of  P^^' 
since  the  subspaces  related  to  the  projections  of  that  method  include  those  of  pO. 

We  partition  an  arbitrary  u/,  G  V**  as 

k 


t=i 


r(3)  :»  .^, «  „f  d(3) 


where  V^   '  is  the  range  of  P-  ' .  The  formulas  for  the  Uh,i  are  given  by 


u 


on      Sli  \  Clz 
"h,i  =  ^  ^1  —  harmonic     on      ^2  \  ^3 

on      CI3 


and  for  2  <  i  <  fc  -  2,  and 


Uh  -  UH,i-i  on      Qi  \  ili+i 

Uh,i  =  ^  hi  —  harmonic     on      f2t+i  \  ^i+2 
0  on      ^i+2 
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Since  the  construction  of  Uh,;,_i  and  u^^h  is  stredghtforwaxd,  we  do  not  provide  details.  We 
note  that  on  each  region  fij  \  fi^+i  only  two  of  the  tenns  differ  from  zero  and  that  u/,  j  =  Uh 

on  d(li+i.  By  applying  the  projection  P,-  '  to  Uh,i,  we  find  that  ti^.j  6  V^'.  The  proof  is 
completed  by  showing  that 

aiuh.i,uh,i)  <  const,  an.\n._^^{uh,Uh). 

This  follows  from  a  variant  of  Lemma  3  and  elementary  considerations.  The  proof  of  the 
lower  botind  is  completed  by  summing  over  t. 
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